Lammska cod prediction grid

Intro

Make an evenly spaced UTM prediction grid with all spatially varying covariates for the diet and the biomass data

pkgs <- c("tidyverse", "tidylog", "sp", "raster", "devtools", "RCurl", "sdmTMB", "terra", "ncdf4", "chron") 

if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
}

invisible(lapply(pkgs, library, character.only = T))

# Source code for map plots
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")

# Packages not on CRAN
# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path
home <- here::here()

Read data and depth-raster

Read stomach data

d <- read_csv(paste0(home, "/data/clean/stomachs.csv")) |>
  mutate(depth_sc = (depth - mean(depth))/sd(depth),
         herring_sc = (herring - mean(herring))/sd(herring),
         saduria_sc = (saduria - mean(saduria))/sd(saduria),          
         sprat_sc = (sprat - mean(sprat))/sd(sprat),
         other_invert_sc = (other_invert - mean(other_invert))/sd(other_invert),
         other_sc = (other - mean(other))/sd(other),
         other_fish_sc = (other_fish - mean(other_fish))/sd(other_fish),
         benth_fish_sc = (benth_fish - mean(benth_fish))/sd(benth_fish),
         year_f = as.factor(year),
         month_f = as.factor(month),
         ices_rect = as.factor(ices_rect),
         pred_length_sc = (pred_length - mean(pred_length)) / sd(pred_length)
) 

glimpse(df)

function (x, df1, df2, ncp, log = FALSE)  

Make the grid with depth

First make a grid for the biomass data, then subset that based on the extend of the stomach data

x <- d$X
y <- d$Y
z <- chull(x, y)

coords <- cbind(x[z], y[z])

coords <- rbind(coords, coords[1, ]) # close the loop

plot(coords[, 1] ~ coords[, 2]) # plot data

sp_poly <- sp::SpatialPolygons(
  list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
  )

sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
                                           data = data.frame(ID = 1)
                                           )
cell_width <- 3 #km?

pred_grid_tmp <- expand.grid(
  X = seq(min(d$X), max(d$X), cell_width),
  Y = seq(min(d$Y), max(d$Y), cell_width),
  year = c(1963:2022)
  )

ggplot(pred_grid_tmp |> filter(year == 2019), aes(X, Y)) +
  geom_point(size = 0.1) +
  theme_void() +
  coord_sf()

filter: removed 3,361,230 rows (98%), 56,970 rows remaining

sp::coordinates(pred_grid_tmp) <- c("X", "Y")

inside <- !is.na(sp::over(pred_grid_tmp, as(sp_poly_df, "SpatialPolygons")))

pred_grid_tmp <- pred_grid_tmp[inside, ]

pred_grid_tmp <- as.data.frame(pred_grid_tmp)

ggplot(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000)) + 
  geom_point(size = 0.001, alpha = 0.5) +
  NULL

filter: removed 2,492,219 rows (98%), 42,241 rows remaining

plot_map +
  geom_point(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
  NULL

filter: removed 2,492,219 rows (98%), 42,241 rows remaining

Warning: Removed 16168 rows containing missing values (`geom_point()`).

# Add lat and lon
# Need to go from UTM to lat long for this one...
# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
xy <- as.matrix(pred_grid_tmp |> dplyr::select(X, Y) |> mutate(X = X*1000, Y = Y*1000))

mutate: changed 2,534,460 values (100%) of 'X' (0 new NA)

        changed 2,534,460 values (100%) of 'Y' (0 new NA)

v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84  +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]

pred_grid_tmp$lon <- lonlat[, 1]
pred_grid_tmp$lat <- lonlat[, 2]

ggplot(filter(pred_grid_tmp, year == 1999), aes(lon, lat)) + geom_point()

filter: removed 2,492,219 rows (98%), 42,241 rows remaining

# Add depth now to remove islands and remaining land
# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
# https://emodnet.ec.europa.eu/geoviewer/
dep_raster <- terra::rast(paste0(home, "/data/Mean depth natural colour (with land).nc"))
class(dep_raster)

[1] "SpatRaster"
attr(,"package")
[1] "terra"

crs(dep_raster, proj = TRUE)

[1] "+proj=longlat +datum=WGS84 +no_defs"

plot(dep_raster)

pred_grid_tmp$depth <- terra::extract(dep_raster, pred_grid_tmp |> dplyr::select(lon, lat))$elevation

ggplot(pred_grid_tmp, aes(lon, lat, color = depth*-1)) + 
  geom_point()

pred_grid_tmp$depth <- pred_grid_tmp$depth*-1 # To make depth from negative elevation in relation to sea surface??

pred_grid_tmp <- pred_grid_tmp |> drop_na(depth)

drop_na: removed 1,014,300 rows (40%), 1,520,160 rows remaining

pred_grid_tmp |> 
  filter(year == 1999) |> 
  drop_na(depth) |> 
  #mutate(water = ifelse(depth < 0.00000001, "N", "Y")) |> 
  ggplot(aes(X*1000, Y*1000, fill = depth)) + 
  geom_raster() +
  NULL

filter: removed 1,494,824 rows (98%), 25,336 rows remaining

drop_na: no rows removed

plot_map + 
  geom_point(data = pred_grid_tmp, aes(X*1000, Y*1000), size = 0.001) + 
  geom_sf() #Simple feature (https://r-spatial.github.io/sf/articles/sf1.html)

Warning: Removed 534240 rows containing missing values (`geom_point()`).

plot_map + 
  geom_raster(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000, fill = depth), size = 0.001) + 
  geom_sf()

filter: removed 1,494,824 rows (98%), 25,336 rows remaining

Warning in geom_raster(data = filter(pred_grid_tmp, year == 1999), aes(X * :
Ignoring unknown parameters: `size`

Warning: Removed 9074 rows containing missing values (`geom_raster()`).

Add ICES areas

VT: changed dat_full to pred_grid_tmp

# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile(paste0(home, "/data/ICES_StatRec_mapto_ICES_Areas/StatRec_map_Areas_Full_20170124.shp"))
head(shape)

  OBJECTID ID ICESNAME SOUTH WEST NORTH EAST AREA_KM2 stat_x stat_y Area_27
1        1 47     47A0  59.0  -44  59.5  -43     3178  -43.5  59.25  14.b.2
2        2 48     48A0  59.5  -44  60.0  -43     3132  -43.5  59.75  14.b.2
3        3 49     49A0  60.0  -44  60.5  -43     3085  -43.5  60.25  14.b.2
4        4 50     50A0  60.5  -44  61.0  -43     3038  -43.5  60.75  14.b.2
5        5 51     51A0  61.0  -44  61.5  -43     2991  -43.5  61.25  14.b.2
6        6 52     52A0  61.5  -44  62.0  -43     2944  -43.5  61.75  14.b.2
          Perc      MaxPer RNDMaxPer AreasList Shape_Leng Shape_Area
1 100.00000000 100.0000000       100    14.b.2          3        0.5
2  84.12674145  84.1267414        84    14.b.2          3        0.5
3  24.99803694  24.9980369        25    14.b.2          3        0.5
4  11.97744244  11.9774424        12    14.b.2          3        0.5
5   3.89717625   3.8971762         4    14.b.2          3        0.5
6   0.28578428   0.2857843         0    14.b.2          3        0.5

pts <- SpatialPoints(cbind(pred_grid_tmp$lon, pred_grid_tmp$lat), 
                     proj4string = CRS(proj4string(shape)))

pred_grid_tmp$subdiv <- over(pts, shape)$Area_27

# Rename subdivisions to the more common names and do some more filtering (by sub div and area)
sort(unique(pred_grid_tmp$subdiv))

 [1] "3.a.20"   "3.a.21"   "3.b.23"   "3.d.24"   "3.d.25"   "3.d.26"  
 [7] "3.d.27"   "3.d.28.1" "3.d.28.2" "3.d.29"   "3.d.32"  

pred_grid_tmp <- pred_grid_tmp |> 
  mutate(sub_div = factor(subdiv),
         sub_div = fct_recode(subdiv,
                              "24" = "3.d.24",
                              "25" = "3.d.25",
                              "26" = "3.d.26",
                              "27" = "3.d.27",
                              "28" = "3.d.28.1",
                              "28" = "3.d.28.2",
                              "29" = "3.d.29"),
         sub_div = as.character(sub_div)) |> 
  filter(sub_div %in% c("24", "25", "26", "27", "28", 2)) |> 
  filter(lat > 54 & lat < 59 & lon < 22)

mutate: new variable 'sub_div' (character) with 10 unique values and 0% NA

filter: removed 376,740 rows (25%), 1,143,420 rows remaining

filter: removed 92,820 rows (8%), 1,050,600 rows remaining

# Add ICES rectangles
pred_grid_tmp$ices_rect <- mapplots::ices.rect2(lon = pred_grid_tmp$lon, lat = pred_grid_tmp$lat)

plot_map +
  geom_raster(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000, fill = depth)) +
  facet_wrap(~sub_div)

filter: removed 1,033,090 rows (98%), 17,510 rows remaining

Warning: Removed 2100 rows containing missing values (`geom_raster()`).

pred_grid <- pred_grid_tmp |> dplyr::select(-subdiv)

Save

# Split in two and save
plot_map +
  geom_raster(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000, fill = depth))

filter: removed 1,033,090 rows (98%), 17,510 rows remaining

Warning: Removed 2100 rows containing missing values (`geom_raster()`).

pred_grid_63_93 <- pred_grid |> filter(year < 1994)

filter: removed 507,790 rows (48%), 542,810 rows remaining

pred_grid_94_22 <- pred_grid |> filter(year > 1993)

filter: removed 542,810 rows (52%), 507,790 rows remaining

saveRDS(pred_grid_63_93, file = paste0(home, "/data/clean/pred_grid_(1_2).rds"))
saveRDS(pred_grid_94_22, file = paste0(home, "/data/clean/pred_grid_(2_2).rds"))

# pred_grid_t <- bind_rows(readRDS(paste0(home, "/data/clean/pred_grid_(1_2).rds")),
#                          readRDS(paste0(home, "/data/clean/pred_grid_(2_2).rds")))
# 
# plot_map +
#   geom_raster(data = filter(pred_grid_t, year == 1999), aes(X*1000, Y*1000, fill = depth))
# 
# 
# saveRDS(pred_grid_63_92, file = paste0(home, "/data/clean/pred_grid_(1_2).rds"))
# saveRDS(pred_grid_93_22, file = paste0(home, "/data/clean/pred_grid_(2_2).rds"))
# 
# # write_csv(pred_grid_63_92, file = paste0(home, "/data/clean/pred_grid_(1_2).csv"))
# # write_csv(pred_grid_93_22, file = paste0(home, "/data/clean/pred_grid_(2_2).csv"))